#include "Utility/Algebra/AlgebraModule.h"

#ifdef SG_USE_MPI

#include "ParAssembly.h"

#include "Utility/Logging/SGLogger.h"

#include "DataStructure/Input/Solver/DOFDirections.h"

#include "DBManager/DBAPI.h"
#include "DBManager/IService/Input/Constraint/ISPCService.h"
#include "DBManager/IService/Input/Load/ISPCDService.h"
#include "DBManager/IService/Input/Node/INodeService.h"
#include "DBManager/IService/Input/Subcase/ISubcaseService.h"

using SG::Algebra::Matrix;
using SG::Algebra::Vector3D;
using SG::DataStructure::Common::Id;
using SG::DataStructure::Common::Real;
using SG::DataStructure::Common::Status;
using SG::DataStructure::FEM::SOLType;
using SG::DataStructure::FEM::DOFDirections;
using SG::DataStructure::FEM::DofSet;
using SG::DataStructure::FEM::Field_Types;
using SG::DataStructure::FEM::NodeData;
using SG::DataStructure::FEM::Subcase;

using namespace SG::DBManager;
using namespace SG::FEMSolver;

/** @brief 调用eigen库实现矩阵分块
 *
 */
void ParAssembly::Matrixchunks (
    const SpMat_t& mat, std::size_t splitIndex, SpMat_t& _OUT mat1, SpMat_t& _OUT mat2, SpMat_t& _OUT mat3, SpMat_t& _OUT mat4)
{
    auto                        M = mat.rows ();
    auto                        N = mat.cols ();

    SG::Algebra::StrideIndexSet row_is1;
    row_is1.first = splitIndex;
    row_is1.n     = M - splitIndex;
    row_is1.step  = 1;

    SG::Algebra::StrideIndexSet row_is2;
    row_is2.first = splitIndex;
    row_is2.n     = M - splitIndex;
    row_is2.step  = 1;

    SG::Algebra::StrideIndexSet col_is1;
    col_is1.first = 0;
    col_is1.n     = 0;
    col_is1.step  = 1;

    SG::Algebra::StrideIndexSet col_is2;
    col_is2.first = splitIndex;
    col_is2.n     = N - splitIndex;
    col_is2.step  = 1;

    // mat1 = mat(0:splitIndex, 0:splitIndex)
    mat.create_submatrix_nosort (mat1, row_is1, col_is1);

    // mat2 = mat(0:splitIndex, splitIndx:N)
    mat.create_submatrix_nosort (mat2, row_is1, col_is2);

    // mat3 = mat(splitIndex:M,0:splitIndex)
    mat.create_submatrix_nosort (mat3, row_is2, col_is1);

    // mat4 = mat(splitIndex:M,splitIndex:N)
    mat.create_submatrix_nosort (mat3, row_is2, col_is2);
}

/** @brief 计算mpc导致的约束方程， um = G *u_n */
void ParAssembly::ComputeMPCControlMatrix ()
{
    SpMat_t Rm_inv;
    m_RmTri.get_inverse (Rm_inv);

    Rm_inv.matrix_matrix_mult (m_RnTri, m_Gm);
    m_Gm.scale (-1.0);
}

ParAssembly& ParAssembly::GetInstance ()
{
    static ParAssembly Instance;
    return Instance;
}

/** @brief 初始化, 使用网格中的所有节点Id, 生成id到节点在数据库索引的映射关系
 *
 * 本函数对节点id进行了升序排列，因此后续查找节点索引可以采用二分查找。</p>
 * 如果要传入新的网格进行计算， 需要调用本类的 `Clear()`接口
 */
void ParAssembly::Initialize (const SG::DBManager::DBServiceFactorySP        dbServiceFactory,
                                const SG::DataStructure::FEM::Subcase&         subcase,
                                const std::set<SG::DataStructure::Common::Id>& loadCases)
{
    // 避免重复初始化
    if (m_isInitialized)
    {
        return;
    }

    // 清空临时变量
    Clear ();
    /** 将 m_isInitialized 显式修改为true。
     * 只有调用 Clear()方法后，再次调用本方法Initialize (), 才会调用后边的语句
     */
    m_isInitialized = true;

    // 计算对用工况的自由度排序
    m_fieldType = Field_Types::FIELD_Elastic;
    if (SOLType::SOL_NLSCSH == subcase.GetSOLType ())
    {
        m_fieldType = Field_Types::FIELD_Thermal;
    }

    const auto nodeIds = DBAPI::GetAllNodeId (dbServiceFactory.get ());
    m_DOFHandler       = SG::FEMSolver::DOFHandler (nodeIds, m_fieldType);

    m_DOFHandler.Compute (dbServiceFactory, subcase);

    // 方程自由度，包括独立自由度和spc约束自由度
    const auto nSparseMatrix = this->GetAcctiveDOFSize ();
    // std::cout << "$$$$$$ nSparseMatrix " << nSparseMatrix << std::endl;
    Vec_t Ff (nSparseMatrix);
    Ff.setZero ();

    const auto nGlobalMSets = m_DOFHandler.GetSlaveDOFSize ();
    Vec_t      Fm(nGlobalMSets);
    Fm.setZero ();

    std::size_t sIndex = 0;
    for (auto loadcase : loadCases)
    {
        m_loadmap[loadcase] = sIndex;
        m_FnSet.push_back (Ff);  // Todo 确认是否是深copy
        m_FmSet.push_back (Fm);
        ++sIndex;
    }
}

Status ParAssembly::FindNodeIndex (SG::DataStructure::Common::Id nodeId, _OUT std::size_t& nodeIndex) const
{
    nodeIndex = m_DOFHandler.FindNodeIndex (nodeId);

    if (-1UL == nodeIndex)
    {
        return Status::STATUS_Failure;
    }

    return Status::STATUS_Success;
}

void ParAssembly::FindNodeIndex (const SG::DataStructure::Common::Id* nodeIds, _OUT std::vector<std::size_t>& nodeIndexes) const
{
    for (auto& index : nodeIndexes)
    {
        index = m_DOFHandler.FindNodeIndex (*nodeIds);
        ++nodeIds;  // 移动指针，指向下一个要处理的节点
    }
}

void ParAssembly::GetGlobalDOFId (const std::vector<std::size_t>&       nodeIndexes,
                                    const SG::DataStructure::FEM::DofSet& DOFs,
                                    _OUT std::vector<int>& dofGlobalIds) const
{
    m_DOFHandler.GetGlobalIndex (nodeIndexes, DOFs, dofGlobalIds);
}

void ParAssembly::GetGlobalDOFId (const std::size_t                            nodeIndex,
                                    const SG::DataStructure::FEM::DOFDirections& dof,
                                    _OUT int&                                    dofGlobalId) const
{
    dofGlobalId = m_DOFHandler.GetGlobalIndex (nodeIndex, dof);
}

/** @brief 查询节点指定自由度的全局编号
 *
 * 首先通过节点id，找到节点的索引， 再读取该节点对应的自由度全局编号
 */
void ParAssembly::ComputeGlobalDOFId (const SG::DataStructure::Common::Id   nodeId,
                                      const SG::DataStructure::FEM::DofSet& DOFs,
                                      _OUT std::vector<int>& dofGlobalIds)
{
    auto nodeIndex = m_DOFHandler.FindNodeIndex (nodeId);
    if (-1 == nodeIndex)
    {
        // 找不到编号为nodeId的节点，将全局自由度全部设为非激活自由度
        for (auto& id : dofGlobalIds)
        {
            id = -1;
        }
        return;
    }

    m_DOFHandler.GetGlobalIndex ({ static_cast<std::size_t> (nodeIndex) }, DOFs, dofGlobalIds);
}

/** @brief 矩阵组装通用过程， 整体刚度矩阵，应力刚度，（质量矩阵）等按照是否做为mpc从自由度组装到不同的矩阵块
 *
 */
auto Assemby = [] (int                     nSparse,
                   int                     nEffectiveDof,
                   const std::vector<int>& globalIndex,
                   const Matrix&           stiffness,
                   SpMat_t&                Kff,
                   SpMat_t&                Kfm,
                   SpMat_t&                Kmm) -> void {
    std::size_t iIndex{ 0 }, jIndex{ 0 };
    for (const auto& iDof : globalIndex)
    {
        // iDof < 0, 该自由度没有激活
        assert (iDof >= 0 && iDof < nEffectiveDof);
        if (iDof < nSparse)
        {
            jIndex = 0;
            for (const auto& jDof : globalIndex)
            {
                assert (jDof > -1 && jDof < nEffectiveDof);
                assert (!std::isnan (stiffness (iIndex, jIndex)));
                if (jDof < nSparse)
                {
                    Kff.add (iDof, jDof, stiffness (iIndex, jIndex));
                }
                else
                {
                    Kfm.add (iDof, jDof - nSparse, stiffness (iIndex, jIndex));
                }
                ++jIndex;
            }
        }
        else
        {
            jIndex = 0;
            for (const auto& jDof : globalIndex)
            {
                assert (jDof > -1 && jDof < nEffectiveDof);
                assert (!std::isnan (stiffness (iIndex, jIndex)));
                // jIndex < nSparse, Kmf不需要组装
                if (jDof >= nSparse)
                {
                    Kmm.add (iDof - nSparse, jDof - nSparse, stiffness (iIndex, jIndex));
                }
                ++jIndex;
            }
        }
        ++iIndex;
    }

    // Ensure that the values are consistent across processors.
    Kff.close ();
    Kfm.close ();
    Kmm.close ();
};

void ParAssembly::AccumulateStiffness (const std::vector<int>& globalIndex, const Matrix& stiffness)
{
    const int nSparse       = static_cast<int> (GetAcctiveDOFSize ());
    const int nEffectiveDof = static_cast<int> (GetAcctiveDOFSize () + m_DOFHandler.GetSlaveDOFSize ());

    Assemby (nSparse, nEffectiveDof, globalIndex, stiffness, m_TriKff, m_TriKfm, m_TriKmm);
}

void ParAssembly::AccumulateCoupledMass (const std::vector<int>& globalIndex, const Matrix& mass)
{
    const int nSparse       = static_cast<int> (GetAcctiveDOFSize ());
    const int nEffectiveDof = static_cast<int> (GetAcctiveDOFSize () + m_DOFHandler.GetSlaveDOFSize ());

    Assemby (nSparse, nEffectiveDof, globalIndex, mass, m_TriCoupledMff, m_TriCoupledMfm, m_TriCoupledMmm);
}

/** @brief 累加应力刚度的贡献 */
void ParAssembly::AccumulateStressStiffness (const std::vector<int>& globalIndex, const Matrix& kSigma)
{
    const int nSparse       = static_cast<int> (GetAcctiveDOFSize ());
    const int nEffectiveDof = static_cast<int> (GetAcctiveDOFSize () + m_DOFHandler.GetSlaveDOFSize ());
    Assemby (nSparse, nEffectiveDof, globalIndex, kSigma, m_TriK4ff, m_TriK4fm, m_TriK4mm);
}

void ParAssembly::AccumulateLoad (SG::DataStructure::Common::Id sid, const std::vector<int>& globalDof, const Vector3D& load)
{
    assert (3 == globalDof.size ());
    const int nSparse       = static_cast<int> (GetAcctiveDOFSize ());
    const int nEffectiveDof = static_cast<int> (GetAcctiveDOFSize () + m_DOFHandler.GetSlaveDOFSize ());
    assert (0 != m_loadmap.count (sid));

    auto& Ff = m_FnSet[m_loadmap[sid]];  // 将变量通过引用绑定
    auto& Fm = m_FmSet[m_loadmap[sid]];

    // 节点对应自由度可能发生越界。（对应自由度可能是约束自由度）

    // "载荷不能加到非激活自由度 && 全局自由度编号不能发生越界
    assert (globalDof[0] >= 0 && globalDof[0] < nEffectiveDof);
    if (globalDof[0] < nSparse)
    {
        Ff.add (globalDof[0], load.m_x);
    }
    else
    {
        Fm.add (globalDof[0] - nSparse, load.m_x);
    }

    assert (globalDof[1] >= 0 && globalDof[1] < nEffectiveDof);
    if (globalDof[1] < nSparse)
    {
        Ff.add (globalDof[1], load.m_y);
    }
    else
    {
        Fm.add (globalDof[1] - nSparse, load.m_y);
    }

    assert (globalDof[2] >= 0 && globalDof[2] < nEffectiveDof);
    if (globalDof[2] < nSparse)
    {
        Ff.add (globalDof[2], load.m_z);
    }
    else
    {
        Fm.add (globalDof[2] - nSparse, load.m_z);
    }

    // Ensure that the values are consistent across processors.
    Ff.close ();
    Fm.close ();
}

void ParAssembly::AccumulateLoad (SG::DataStructure::Common::Id sid, const std::vector<int>& globalDof, const Matrix& load)
{
    assert (load.GetRowSize () == globalDof.size ());
    assert (1 == load.GetColSize ());
    assert (0 != m_loadmap.count (sid));

    auto& Ff = m_FnSet[m_loadmap[sid]];  // 将变量通过引用绑定
    auto& Fm = m_FmSet[m_loadmap[sid]];

    const int nSparse       = static_cast<int> (GetAcctiveDOFSize ());
    const int nEffectiveDof = static_cast<int> (GetAcctiveDOFSize () + m_DOFHandler.GetSlaveDOFSize ());

    auto nComponent = load.GetRowSize ();
    for (std::size_t index{ 0 }; index < nComponent; ++index)
    {
        const auto globalDofId = globalDof[index];
        const auto fi          = load (index, 0);
        // "载荷不能加到非激活自由度 && 全局自由度编号不能发生越界
        if (globalDofId < 0 || globalDofId >= nEffectiveDof)
        {
            if (fabs (fi) > 1.0e-12)
            {
                SGWarn << __FILE__ << ":" << __LINE__ << "invaid load component";
            }
            continue;
        }

        // 按照自由度是否为 mpc从自由度分别组装
        if (globalDofId < nSparse)
        {
            Ff.add (globalDofId, fi);
        }
        else
        {
            Fm.add (globalDofId - nSparse, fi);
        }
    }

    // Ensure that the values are consistent across processors.
    Ff.close ();
    Fm.close ();
}

/** 累加约束方程
 * @param [in] gdofn 方程系数列 对应的全局自由度
 * @param [in] gdofm 方程系数行 对应的全局自由度， 约束方程对哪一个约束方程又贡献
 *
 */
void ParAssembly::AccumulationControlMatrix (const std::vector<int>& gdofn, const std::vector<int>& gdofm, const Matrix& RnRm)
{
    assert (gdofn.size () == RnRm.GetColSize ());
    assert (gdofm.size () == RnRm.GetRowSize ());

    const auto nActiveDofs  = static_cast<int> (GetAcctiveDOFSize ());             // 网格中自由的自由度数量
    const auto nGlobalMSets = static_cast<int> (m_DOFHandler.GetSlaveDOFSize ());  // 网格中多点约束方程的数量
    const auto nDofM        = gdofm.size ();
    for (std::size_t row = 0; row < nDofM; ++row)
    {
        // 对应元素组装到约束方程的行号, 即组装到哪一个约束方程中
        auto iEquation = gdofm.at (row) - nActiveDofs;  //
        assert (iEquation >= 0);

        // 依次处理矩阵的列
        std::size_t col{ 0 };
        for (auto dof : gdofn)
        {
            if (dof < 0)  // 主自由度为非激活自由度
            {
                SGError << " 主自由度必须为为自由状态" << std::endl;
                continue;
            }
            assert (dof < (nActiveDofs + nGlobalMSets));

            const auto& val = RnRm (row, col);
            assert (!std::isnan (val));

            // 避免将 0.0填入控制方程，导致影响控制方程的稀疏性
            // 测试发现此次过滤0.0的实现存在问题
            // if (fabs (val) > 1.0e-18)
            // {
            if (dof < nActiveDofs)
            {
                m_RnTri.add (iEquation, dof, val);
            }
            else
            {
                m_RmTri.add (iEquation, dof - nActiveDofs, val);
            }
            // }

            ++col;
        }
    }

    // Ensure that the values are consistent across processors.
    m_RnTri.close ();
    m_RmTri.close ();
}

void ParAssembly::ClearBCLoad ()
{
    m_BCLoadIdSet.clear ();

    m_FnBCSet.clear ();
    m_FmBCSet.clear ();
    m_FfBCSet.clear ();
}

void ParAssembly::AccumulateBCLoad (SG::DataStructure::Common::Id sid, const std::vector<int>& globalDof, const SG::Algebra::Matrix& load)
{
    assert (load.GetRowSize () == globalDof.size ());
    assert (1 == load.GetColSize ());
    if (0 == m_BCLoadIdSet.count (sid))
    {
        // 激活当前工况
        m_BCLoadIdSet.insert (sid);
        // 方程自由度，包括独立自由度和spc约束自由度
        const auto nSparseMatrix = this->GetAcctiveDOFSize ();
        // std::cout << "$$$$$$ nSparseMatrix " << nSparseMatrix << std::endl;
        Vec_t Ff (nSparseMatrix);
        Ff.setZero ();
        m_FnBCSet[sid]          = Ff;
        const auto nGlobalMSets = m_DOFHandler.GetSlaveDOFSize ();
        Vec_t      Fm (nGlobalMSets);
        Fm.setZero ();
        m_FmBCSet[sid] = Ff;
    }

    auto& Ff = m_FnBCSet[sid];  // 将变量通过引用绑定
    auto& Fm = m_FmBCSet[sid];

    const int nSparse       = static_cast<int> (GetAcctiveDOFSize ());
    const int nEffectiveDof = static_cast<int> (GetAcctiveDOFSize () + m_DOFHandler.GetSlaveDOFSize ());

    auto nComponent = load.GetRowSize ();
    for (std::size_t index{ 0 }; index < nComponent; ++index)
    {
        const auto globalDofId = globalDof[index];
        const auto fi          = load (index, 0);
        // "载荷不能加到非激活自由度 && 全局自由度编号不能发生越界
        if (globalDofId < 0 || globalDofId >= nEffectiveDof)
        {
            if (fabs (fi) > 1.0e-12)
            {
                SGWarn << __FILE__ << ":" << __LINE__ << "invaid load component";
            }
            continue;
        }

        // 按照自由度是否为 mpc从自由度分别组装
        if (globalDofId < nSparse)
        {
            Ff.add (globalDofId, fi);
        }
        else
        {
            Fm.add (globalDofId - nSparse, fi);
        }
    }

    // Ensure that the values are consistent across processors.
    Ff.close ();
    Fm.close ();
}

void applyAutoSPC (const DBServiceFactorySP dbServiceFactory,
                   const std::vector<int>&  dofGlobalIndexes,
                   const int                kDim,
                   const std::size_t        nodeIndex,
                   const Id                 nodeId,
                   Real                     maxValue,
                   SpMat_t& _OUT            K,
                   bool                     isRot)
{
    // @todo: complement the required implementations
    // ...
}

/** @brief 读取并组装 spc对应的位移矩阵
 *
 * 对于GRID的“PS”, 由于位移值为0.0. 因此不需要处理有关的值
 * @param      dbServiceFactory 数据库
 * @param      spcId     单点约束集合编号
 * @param      nTotalSet 模型中所有激活自由度数量
 * @param      nfSet     独立自由度数量
 * @param[out] Ys        spc对应的位移矩阵, 向量应该提取指定相应的size, 并且将所有元素初始化为0.0;
 */
void assemblyYs (const SG::DBManager::DBServiceFactorySP dbServiceFactory,
                 const SG::FEMSolver::DOFHandler&        DOFHandler,
                 Subcase                                 subcase,
                 int                                     nTotalSet,
                 int                                     nfSet,
                 _OUT Vec_t&                       Ys)
{
    const auto nSDOFSet = nTotalSet - nfSet;
    Ys.resize (nSDOFSet);
    Ys.setZero ();

    /** spc 和spc1 数据 */
    SG::DataStructure::FEM::SPCDataArray spcSet;
    dbServiceFactory->get<ISPCService> ()->FindBySId (subcase.GetSPCId (), _OUT spcSet);
    /** 叠加spcd数据 */
    SG::DataStructure::FEM::SPCDataArray spcdSet;
    dbServiceFactory->get<ISPCDService> ()->FindBySId (subcase.GetLoadId (), spcdSet);
    spcSet.insert (spcSet.end (), spcdSet.begin (), spcdSet.end ());

    for (auto& spc : spcSet)
    {
        auto nodeIndex = DOFHandler.FindNodeIndex (spc.m_g);
        if (-1 == nodeIndex)
        {
            std::cout << "warning, shouldn't constraint this node [ " << spc.m_g << " ]" << std::endl;
            continue;
        }

        DofSet     indexes = spc.GetIndexes ();
        const auto dofSize = indexes.size ();
        for (std::size_t index = 0; index < dofSize; ++index)
        {
            const auto posi = DOFHandler.GetGlobalIndex (nodeIndex, indexes[index]);
            // 跳过非激活自由度
            if (posi < 0)
            {
                std::cout << "shouldnot constraint " << spc.m_g << " direction " << indexes[index] << std::endl;
                continue;
            }

            // 仅仅需要处理非零情况
            if (fabs (spc.m_d) > 1.0e-10)
            {
                Ys.set (posi - nfSet, spc.m_d);
            }
        }
    }

    // Ensure that the values are consistent across processors.
    Ys.close ();
}

/** @brief 约束处理， 消去mpc的非独立自由度， 处理spc约束，添加autospc计算功能
 *
 * 使用线性消去处理mpc方程， 使用乘大数法处理spc计算功能， 自动进行autoSPC 处理
 */
void ParAssembly::Assembly (const DBServiceFactorySP dbServiceFactory, const SG::DataStructure::Common::Id subcaseId)
{
    SGDebug << "== assembly start ";

    auto       n     = GetAcctiveDOFSize ();
    auto       nDofM = m_DOFHandler.GetSlaveDOFSize ();
    const auto nKDim = static_cast<int> (m_DOFHandler.ComputeSparseMatrixSize ());

    // 找到对角元中的绝对值最大值
    Real  maxValue{ 0.0 };
    Vec_t diag;
    m_K.get_diagonal (diag);
    diag.linfty_norm (maxValue);

    // 处理多点约束方程
    if (NeedComputeControlMatrix ())
    {
        SGDebug << "== deal with Milti Point Constraint ";
        // 计算 m_Gm
        ComputeMPCControlMatrix ();
        SpMat_t GmT;
        m_Gm.get_transpose (GmT);

        SpMat_t Knm_T;
        m_TriKfm.get_transpose (Knm_T);

        SpMat_t mat1;
        GmT.matrix_matrix_mult (Knm_T, mat1);

        SpMat_t mat2;
        mat1.get_transpose (mat2);

        SpMat_t mat3;
        m_Gm.matrix_matrix_mult (m_TriKmm, mat3);

        SpMat_t mat4;
        mat3.matrix_matrix_mult (m_Gm, mat4);

        m_K = m_K + mat1 + mat2 + mat4;

        // 质量矩阵多点约束处理
        if (m_needMass)
        {
            SpMat_t mat5;
            m_TriCoupledMfm.get_transpose (mat5);

            SpMat_t mat6;
            GmT.matrix_matrix_mult (mat5, mat6);

            SpMat_t mat7;
            m_TriCoupledMmm.matrix_matrix_mult (m_Gm, mat7);

            SpMat_t mat8;
            GmT.matrix_matrix_mult (m_TriCoupledMmm, mat8);

            SpMat_t mat9;
            mat8.matrix_matrix_mult (m_Gm, mat9);

            m_CoupledM = m_CoupledM + mat6 + mat9;
        }

        // 载荷右端项多点约束处理
        const auto nLoadCase = m_FnSet.size ();
        for (std::size_t iLoad{ 0 }; iLoad < nLoadCase; ++iLoad)
        {
            m_FnSet[iLoad] = m_FnSet[iLoad]  + GmT * m_FmSet[iLoad];
        }

        // 右端项特殊边界考虑多点约束
        // const auto nBC = m_BCLoadIdSet.size();
        for (auto sid : m_BCLoadIdSet)
        {
            m_FnBCSet[sid] = m_FnBCSet[sid] + GmT * m_FmBCSet[sid];
        }
    }

    /**  处理单点约束（包括非零强制位移）
     */
    SGDebug << "== deal with Single Point Constraint ";
    auto    para = dbServiceFactory->get<ISubcaseService> ()->FindSubCase (subcaseId);
    Subcase subcase (para);
    assemblyYs (dbServiceFactory, m_DOFHandler, subcase, n, nKDim, _OUT m_Ys);

    /** 将稀疏矩阵按照（1）被spc约束或者（2）做为独立自由度这两种情况 进行分块 */
    SpMat_t     ksf;  // dummy
    SpMat_t     kss;  // dummy
    ParAssembly::Matrixchunks (m_K, nKDim, _OUT m_K, _OUT m_kfs, _OUT ksf, _OUT kss);

    if (m_needMass)
    {
        SpMat_t     msf;  // dummy
        SpMat_t     mss;  // dummy
        SpMat_t     mfs;
        ParAssembly::Matrixchunks (m_CoupledM, nKDim, _OUT m_CoupledM, _OUT mfs, _OUT msf, _OUT mss);
    }

    SG::Algebra::StrideIndexSet idx;
    idx.first = 0;
    idx.step  = 1;
    idx.n     = nKDim;

    /** 右端项约束处理 */
    auto nLoadCase = m_FnSet.size ();
    m_FfSet.resize (nLoadCase);
    for (std::size_t iLoad{ 0 }; iLoad < nLoadCase; ++iLoad)
    {
        auto& Fn = m_FnSet[iLoad];
        // 分块修正单点约束（spc）自由度
        Vec_t h;
        Fn.create_subvector (h, idx);
        m_FfSet[iLoad] = h - m_kfs * m_Ys;
    }

    // 右端项特殊边界考虑多点约束
    // const auto nBC = m_BCLoadIdSet.size();
    for (auto sid : m_BCLoadIdSet)
    {
        auto& Fn       = m_FnBCSet[sid];
        Vec_t h;
        Fn.create_subvector (h, idx);
        m_FfBCSet[sid] = h - m_kfs * m_Ys;
    }

    // auto SPC
    if (Field_Types::FIELD_Elastic == m_fieldType)
    {
        static DofSet    transDOFSet{ DOFDirections::DOF_Ux, DOFDirections::DOF_Uy, DOFDirections::DOF_Uz };
        static DofSet    rotDOFSet{ DOFDirections::DOF_Rotx, DOFDirections::DOF_Roty, DOFDirections::DOF_Rotz };
        auto             nNode = m_DOFHandler.GetNodeSize ();
        std::vector<int> dofGlobalIndexes (3);
        const auto&      nodeIds = m_DOFHandler.GetNodeIds ();
        for (std::size_t iNode = 0; iNode < nNode; ++iNode)
        {
            m_DOFHandler.GetGlobalIndex ({ iNode }, transDOFSet, _OUT dofGlobalIndexes);
            auto nodeId = nodeIds[iNode];
            applyAutoSPC (dbServiceFactory, dofGlobalIndexes, nKDim, iNode, nodeId, maxValue, _OUT m_K, false);
            // rot
            m_DOFHandler.GetGlobalIndex ({ iNode }, rotDOFSet, _OUT dofGlobalIndexes);
            applyAutoSPC (dbServiceFactory, dofGlobalIndexes, nKDim, iNode, nodeId, maxValue, _OUT m_K, true);
        }
    }

    // 压缩矩阵， 去掉 （_， _）占位元素
    m_K.close();
}

void ParAssembly::ComputeBVector (SG::DataStructure::Common::Id loadCase, _OUT Vec_t& b) const
{
    auto it = m_loadmap.find (loadCase);
    if (it == m_loadmap.end ())
    {
        const auto nKDim = static_cast<int> (m_DOFHandler.ComputeSparseMatrixSize ());
        b.resize (nKDim);
        b.setZero ();
    }
    else
    {
        b = m_FfSet[it->second];
    }

    for (auto sid : m_BCLoadIdSet)
    {
        b = b + m_FfBCSet.at (sid);
    }
}

void ParAssembly::ComputeBVector (const SG::DataStructure::FEM::LoadCombination& load, _OUT Vec_t& b) const
{
    const auto nSparseMatrix = m_DOFHandler.ComputeSparseMatrixSize ();
    b.resize (nSparseMatrix);
    b.setZero ();

    for (auto item : load.m_Items)
    {
        auto loadCase = item.m_subSet;
        auto it       = m_loadmap.find (loadCase);
        if (it == m_loadmap.end ())
        {
            assert (0);
            return;
        }
        const auto index = it->second;
        b = b + item.m_subfactor * m_FfSet[index];
    }

    b = load.m_facotor * b;
}

const SpMat_t ParAssembly::ComputeStressStiffnessMatrix ()
{
    auto nNSet = GetAcctiveDOFSize ();
    auto nMSet = m_DOFHandler.GetSlaveDOFSize ();

    SpMat_t K4 (nNSet, nNSet);

    return K4;
}

/**
 *  Real 类型的位移恢复
 */
template <>
Pre_DECLSPEC void ParAssembly::RecoveryDisplacement (Vec_t& displacement, _OUT Vec_t& extendDisplacement) const
{
}

/**
 * std::complex<double> 类型的位移恢复
 */
template <>
Pre_DECLSPEC void ParAssembly::RecoveryDisplacement (ComplexVec_t& displacement, _OUT ComplexVec_t& extendDisplacement) const
{

}

#endif  // SG_USE_MPI